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simple,  hence  the  computer  simulation  of  man-in-the-loop  AAA  tracking 
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is  concluded  that  the  antiaircraft  gunner  model  based  on  the  observer 
theory  can  be  accurately  and  efficiently  used  to  study  AAA  weapon  effective¬ 
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SECTION  I 


INTRODUCTION 

The  Manned- Systems  Effectiveness  Division  of  the  Aerospace  Medical 
Research  Laboratory  (AMRL/ME) ,  WPAFB,  Ohio,  performs  research  on  the  gunner's 
tracking  response  in  Anti-Aircraft  Artillery  (AAA)  systems.  To  analyze  the 
performance  of  an  AAA  system,  a  mathematical  gunner  model  representing  the 
control  characteristics  of  the  gunner  in  the  compensating  tracking  task  is 
required.  Several  human  operator  models  have  been  developed: 

1.  McRuer  Crossover  Model  [1]  -  based  on  classical  control  theory, 

2.  Optimal  Control  Model  [2],  [3],  [4]  -  based  on  optimal  control  and 
estimation  theory, 

3.  PID  Structure  Modified  Optimal  Control  Model  [5]  -  Simplified  Optimal 
Control  Model. 

A  brief  comparison  of  these  models  can  be  found  in  Appendix  A.  Although 
most  of  these  models  can  predict  tracking  errors,  they  either  have  a  complicated 
model  structure  and/or  it  is  difficult  to  determine  values  for  the  model 
parameters  for  a  given  weapon  system. 

We  wish  to  design  a  model  whose  structure  is  simple,  whose  parameters  could 
be  identified  systematically,  whose  computer  implementation  would  be  fast  and 
efficient,  and  whose  output  would  accurately  describe  the  tracker's  response 
characteristics. 

This  report  presents  an  AAA  gunner  model  based  on  the  Luenberger  observer 
theory  [6],  [7],  [8].  The  observer  theory  provides  a  new  method  to  obtain  an 
approximate  estimation  for  the  state  of  an  observed  syscem.  The  characteristics 
of  an  observer  are  somewhat  free  to  the  extent  that  they  can  be  determined  by 
the  designer  through  the  proper  selection  of  an  observer  gain.  There  is  no 
Ricrati  equation  involved  in  the  observer  design.  The  simplicity  of  the  observer 
design  and  its  capability  for  state  estimation  make  the  observer  theory  an 
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attractive  design  method.  Then  the  estimated  state  vector  can  be  used  to 
implement  a  linear  state  variable  feedback  controller  which  represents  the 
gunner's  tracking  function.  Furthermore,  it  is  assumed  that  the  effects  of 
all  randomness  sources  in  the  AAA  closed  loop  system  can  be  lumped  into  one 
remnant  element.  Therefore,  the  structure  of  this  model  is  simple.  In  addition, 
a  parameter  identification  program  based  on  least  squares  curve-fitting  method 
and  the  Gauss  Newton  algorithm  [9]  was  developed  for  this  model.  Hence,  its 
parameters  can  be  easily  determined.  A  computer  simulation  program  OMS 
(Observer  Model  Simulation)  of  the  AAA  tracking  task  using  this  model  was  also 
developed.  Simulation  results  showed  that  the  model  predictions  of  the  tracking 
errors  were  in  excellent  agreement  with  the  actual  gunner  response  data  of  the 
manned  AAA  simulation  conducted  at  the  Aerospace  Medical  Research  Laboratory, 
WPAFB,  Ohio.  The  computer  execution  time  of  the  AAA  simulation  using  this 
simple  model  is  very  fast.  The  description  of  a  AAA  gun  system  and  the  design 
of  the  observer  model  are  included  in  Section  II.  Section  III  describes  the 
method  to  determine  the  model  parameters.  Discretization  techniques  for  computer 
implementation  will  be  included  in  Section  IV  along  with  the  computer  simulation 
results.  The  conclusion  will  be  discussed  in  Section  V. 


SECTION  II 
MODEL  DEVELOPMENT 

A.  Description  of  An  AAA  Gun  System 

The  tracking  task  of  an  anti-aircraft  artillery  (AAA)  gun  system  can  be 
described  by  a  closed  loop  block  diagram  as  shown  in  Figure  1.  Two 
gunners,  one  each  for  azimuth  and  elevation  axes,  play  the  role  of 
controller  in  the  man-machine  feedback  control  system.  From  the  visual 
display,  each  gunner  observed  the  tracking  error  ,  which  is  the  differ¬ 
ence  between  the  target  position  angle  0T,  and  the  gunsight  angle  0  . 
Independently,  the  gunners  operated  hand  cranks  to  control  the  gunsight 
system  to  align  the  gunsight  angle  (output)  with  the  target  position 
angle  (input).  Therefore,  the  azimuth  tracking  task  is  decoupled  from 
the  elevation  tracking  task  in  this  AAA  system.  Four  trajectories,  Figure 
2,  of  the  target  aircraft  were  used  as  input  to  the  AAA  system.  These 

trajectories  are  deterministic  functions  of  time,  although  the  gunners 

•  •  • 

do  not  know  their  dynamic  properties,  0T,  0T>  etc.  In  order  to  develop 
a  mathematical  model  of  the  gunner  response  characteristics  in  an  AAA 
compensatory  tracking  task,  we  first  need  to  describe  the  mathematical 
representations  of  the  gunsight  dynamics.  Let  us  consider  an  AAA  system 
with  the  following  gunsight  transfer  function. 


9  (s) 

8 


64  (s  +  1) 


U  (s)  s(s4  +  12s  +  64) 

Based  on  frequency  domain  analysis  of  the  inputs  (target  trajectories), 
it  was  found  that  the  frequency  bandwidths  of  all  the  trajectories  in 
Figure  2  are  around  0.2  Hz.  Thus,  a  simplified  gunsight  transfer 

function 

9g  (s)  ,  1 

U  (s)  s 

can  be  used  for  the  model  design  and  simulation  analysis. 
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(1) 


FIGURE  1 :  BLOCK  DIAGRAM  OF  AN  AAA  CLOSED  LOOP  SYSTEM 


Equation  (1)  is  equivalent  to  6  (t)  =  u(t).  Furthermore,  Figure  1  defines 

o 

the  tracking  error  eT(t)  to  be  (t >  -  9^(t).  The  tirae  derivative  of  eT(t) 
is: 

e_(t)  -  0  (t)  -  6  (t) 

T  .T  g 

=  eT(t>  -  u(t) 

Now  let  us  introduce  state  variables 

A  • 

xL(t)  =  eT(t) 

x2(t)  =  eT(t)  =  0T(t)  -  o  (t) 

The  time  derivatives  of  these  state  variables  are 

ix(t)  a  0T(t) 

x2(t)  =  Xl(t)  -  u(t) 

T 

Let  x(t)  be  [x^(t),  x2(t)]  .  The  state  space  equation  of  the  gunsight 
dynamics  and  target  motion  can  be  expressed  as 

x(t)  =  Ax(t)  +  Bu (t )  +  F0^(t)  (2) 

where  A,  B,  and  F  are  matrices  defined  as  follows: 


0 

0 

»  B  = 

1 - 

o 

I _ 

A  = 

1 

O 

1 

L-M 

The  scalars  u  and  0T  are  the  control  output  of  an  AAA  gunner  and  the  target 
acceleration  respectively. 

The  tracker's  observation  of  the  tracking  error  is  represented  in  the 
measurement  equation: 

y(t)  =  Cx(t)  (3) 

where  C  is  a  row  vector  [01]. 

B.  Human  Operator's  Internal  Model 

To  develop  the  observer  model  for  the  human  operator's  tracking  response, 
it  is  necessary  to  obtain  his  internal  model  of  the  controlled  plant  and 
target  motion,  Figure  3.  This  internal  model  describes  the  understanding 
or  knowledge  the  gunner  (i.e.  human  operator)  has  about  the  real  system 


V-  — 
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(including  gunsight  dynamics,  display,  control  crank,  and  target  motion). 

For  simplicity,  it  is  assumed  that  the  human  operator's  internal  model  of 
gunsight  dynamics  is  identical  to  the  mathematical  model  of  gunsight 
dynamics. 

As  mentioned,  the  tracker  has  information  about  target  trajectory  6 

1  y  ' 

and  maybe  target  velocity  9T,  but  not  precise  information  of  the  target 
acceleration  0  .  Hence,  the  last  term,  0T>  in  Eq.  (1)  will  not  be  included 
in  the  human  operator's  internal  model  of  the  system.  This  is  the  target 
uncertainty  problem  whose  effect  will  be  included  in  a  remnant  element, 
considered  in  the  next  section. 

In  addition,  the  observation  noise  which  might  be  associated  with 
measured  output  y  will  not  be  considered  explicitly.  Its  random  effect  will 
also  be  included  in  the  remnant  element  for  simplicity.  Therefore,  the 
human  operator's  internal  model  of  the  system  can  be  described  as  follows: 

x  (t)  =  Ax, (t)  +  Bu (t) 

1  1  (4) 

y  (t)  =  C x(t) 

The  state  variable  x^(t)  of  the  internal  model  also  has  dimension  two. 

'  C.  Observer  Model 

An  observer  is  itself  a  dynamic  system  whose  function  is  to  reconstruct  the 
state  variable  of  a  given  system  (in  this  case,  the  human  operator's  internal 
model).  An  identity  observer  is  an  observer  which  has  the  same  dynamic  order 
as  the  observed  system.  In  this  report,  an  identity  observer  is  designed 
to  estimate  the  states  of  the  gunsight  and  the  target  motion.  The  structure  of 
an  identity  observer  used  in  our  observer  model  design  is: 

£(t)  =  Az.(t)  +  Bu^  (t)  +  K  [y(t)  -  Cz.(t)] 

=  (A-KC)z_(t)  +  Buc(t)  +  Ky  (t)  (5) 

z,(t)  is  the  state  variable  of  the  observer  with  dimension  two.  This  vector 
represents  the  estimated  value  of  the  state  variable  x(t)  and  is  used  to  fulfill 


the  state  variable  feedback  controller.  Matrices  A,B,  and  C  are  as  previously 
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described.  K  is  the  observer  gain  vector.  The  dynamic  response  of  the 
observer  is  determined  by  the  matrix  A-KC.  K  is  designed  such  that  the 
eigenvalues  of  A-KC  have  negative  real  parts  (i.e.  A-KC  is  a  stable  matrix). 
Then  the  state  of  the  observer  will  converge  to  the  state  of  the  observed 
system.  According  to  observer  theory,  there  always  exists  a  gain  K  such  that 
A-KC  is  stable  if  the  system  is  observable.  The  definition  of  an  observable 
system  can  be  formed  in  [10].  It  can  be  shown  that  the  system  described  by 
Eq.  (4)  is  observable.  u^(t)  in  Eq. (5)  is  the  feedback  controller  and  is 
defined  as  follows: 

u  (t)  =  -r«(t)  (6) 

c  — 

where  T  is  the  controller  gain,  a  row  vector  with  two  elements.  In  Eqs.  (5) 
and  (6),  the  observer  gain  K  and  the  controller  gain  T  will  be  determined  by 
a  parameter  identification  program. 

The  relation  between  the  control  output  u(t)  of  the  human  operator  and 

u  (t)  is: 
c 

u(t)  =  u  (t)  +  v(t) 
c 

where  v(t)  is  a  random  process  called  remnant  element  which  represents  all 
the  gunner-induced  noises,  (e.g.,  the  effect  of  target  uncertainty,  the 
observation  noise,  the  neuromotor  noise,  etc.)  and  modeling  errors.  The 
idea  to  lump  all  the  random  effects  into  the  remnant  element  is  to  simplify 
the  structure  of  the  model.  The  statistical  properties  of  the  remnant  v(t) 
are : 


and 


E[v(t)]  =  0  for  all  t, 

E [ v ( t )  vT(t)J  =  V(t)6(t-T)  for  all  t  and  T. 


(7) 


where  E [ * ]  denotes  the  expectation  value  of  6  is  the  Dirac  delta  function 
and  V(t)  is  a  function  of  time  to  be  described  in  section  III.  A  block  diagram 
of  the  structure  of  the  observer  model  is  shown  in  Figure  4.  Let  £(t)  be  the 
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FIGURE  4:  BLOCK  DIAGRAM  OF  OBSERVER  MODEL  STRUCTURE 


estimation  error,  i.e.,  £(t)  =  x(t)  ~  2.(t),  then 

e(t)  =  x(t)  -  i(t) 
and  from  Eqs.  (2)  and  (5) 

e(t)  =  (A-KC)e^(t)  +  Bv(t)  +  F9T(t)  (8) 

K  will  be  selected  to  make  A-KC  a  stable  matrix.  Hence,  <e(t)  will  decrease 
exponentially  to  zero  with  time.  The  mathematical  model  of  the  AAA  closed 
loop  system  which  is  composed  of  Eqs.  (2)  and  (5)  can  be  written  as: 

x(t)  =  (A-Br)jc(t)  +  Bre(t)  +  Bv(t)  +  F8T(t) 

(9) 

e(t)  =  (A-KC)e(t)  +  Bv(t)  +  F0^,(t) 


z^t)  =  x(t) 


A,-  A-Br  6r  ,  B, 
0  A-KC 


md  F, 


^(t)  =  A^Ct)  +  BjvCt)  +  F^ft)  (1°) 

Note  that  the  system  matrix  A^  of  the  above  overall  system  is  a  triangular 
matrix.  Hence,  the  eigenvalues  of  this  triangular  matrix  are  the  eigenvalues 
of  matrices  A-BT  and  A-KC.  Now,  if  we  choose  a  proper  control  gain  matrix  V 
and  observer  gain  matrix  K  to  make  A-Br  and  A-KC  stable  matrices  then  the 
overall  system  Is  stable.  Furthermore,  it  can  be  shown  that  the  design  of 
observer  and  the  design  of  controller  can  be  done  separately.  This  simplifies 
the  design  procedures.  The  mean  of  Eq.  (10)  is: 

z^t)  =  AjZ^t)  +  Fn8T(t)  (-q) 

where  z^(t)  =  E[z^(t)]. 

The  covariance  of  this  equation  is: 

P(t)  =  A1P(t)  +  +  B1V(t)B1T  (12) 

where  P(t)=  E[(z^(t)  -  £^(t))(£^(t)  -  z^(t))T].  The  derivation  of  Eq.  (12) 


can  be  found  in  Reference  [11]. 
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SECTION  III 


PARAMETER  IDENTIFICATION  METHOD 

The  structure  of  the  gunner  model  has  been  designed,  To  implement  this 
work,  we  need  to  determine  the  values  of  the  model  parameters.  A  systematic 
procedure  has  been  developed  (Figure  5).  It  involves  ar.  identification  program 
based  on  the  least  squares  curve-fitting  method  and  the  Gauss-Newton  gradient 
algorithm  as  described  in  Appendix  B.  One  requirement  of  this  technique  is  to 
select  a  criterion  function  J(a)  of  the  unknown  parameters  a  to  evaluate  the 
"goodness  of  fit"  between  the  model  predictions  from  Eqs.  01)  and  (12)  and  the 
empirical  data  obtained  from  the  manned  AAA  simulation.  Parameter  identification 
occurs  iteratively  by  converging  to  a  set  of  values  for  which  the  curves  are 
"reasonably"  matched. 

Empirical  data  of  the  tracking  error,  e^,(t),  was  collected  at  the  Aerospace 
Medical  Research  Laboratory,  WPAFB,  Ohio,  from  their  manned  simulator  using 
simulated  target  trajectory  0^  as  input,  and  SDT  represent  the  sample  ensemble 
mean  and  standard  deviation  of  the  tracking  errors  over  sixteen  simulation  runs 
for  a  given  team.  In  Eq.  (11)  the  model  prediction  of  the  ensemble  mean  of  the 
tracking  error,  e^  is  located  in  the  second  component  of  z^(t).  It  can  also  be 
shown  that  the  square  root  of  the  second  diagonal  element  of  the  covariance  matrix 
P(t)  in  Eq.  (12)  is  the  predicted  ensemble  standard  deviation  of  the  tracking 
error.  The  parameter  identification  was  done  in  two  parts  (curve-fitting  of  the 
ensemble  mean  and  ensemble  standard  deviation  of  the  tracking  error).  For  both 
parts,  angular  information  was  input  only  from  Trajectory  A  of  Figure  2. 

In  the  first  part,  e^  and  e^,  are  compared  to  determine  the  observer  gain 
K  =  [k^,  k£]  and  the  controller  gain  T  =  [y^,  Yj]-  This  comparison  was  actually 
computed  in  frequency  domain.  The  power  spectral  density  function  (PSD)  of  the 
sample  ensemble  mean  tracking  error  was  calculated,  e^,(u)).  And  the  PSD,  e^,(o))  » 
of  e^  can  be  computed  from  the  state  representation  of  z^(t)  in  Eq.  (11).  The 
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AAA  SIMULATOR 


solution  to  this  differential  equation  is: 


z  (t)  =  r  *(t-T)F  0T(T)dT 
— 1  -00  1  T 


where  <i>( t)  =  EXP(A  t)  is  the  transition  matrix. 
1 


e'(t)  =  r  [<J>  (t-T)  +  <f>  (t-T )  ]  0  (t  )dx  u<u 

T  -“2  1  2  3  T 

t  h  - 

where  $„(t)  is  the  ij  element  of  the  matrix  <f>(t).  The  PSD  of  e^,(t),  denoted  by 

S_ 

e^,(to)  can  be  shown  [10]  to  be 


Se’  (w)  -  | G  (ju»|  S0T(a))  +  | G  (300) )  S0  <u>) 

I  21  i  23  *- 


2  Q,. 

i  ^  a 


=  [  | G  (jto)  j  +  |G  (ji>)|  ]  \(u) 


S"  "  th 

where  0^,(oO  is  the  PSD  of  0,j,(t)  and  G„  (jco)  is  the  ij  element  of  the  square 

matrix  G ( joo)  which  is  defined: 


G(s)  =  L  (4>(t)  ]  (16) 

i.e.,  G(s)  is  the  Laplace  transform  of  the  transition  matrix  <j>(t)  and  s  is  the 
variable  of  the  complex  plane.  It  can  be  shown  that  G(s)  =  (sI-A^)  1  hence  can 


be  rewritten  as: 


C(s)  = 


adj  (sI-A  ) 
_ 1 _ 

det  (sI-A  ) 
l 


where  I  is  the  identity  matrix,  det  and  adj  denote  the  determinant  and  adjoint 
of  a  matrix  respectively. 

Furthermore,  because  of  the  structure  of  A  ,  it  can  be  shown  that  (see 


reference  [8]) 


G(s)  = 


adj  (sI-A^ ) 

det  [sI-(A-BO)  •  det  [sI-(A-KC)] 


Hence 


G  (s)  =  1  +  Yi 

21  s(s-y  ) 

2 

y  s  +  y  k  +  Y 

g23(s)  =  1  ]  2  2 

(s-Y  )  (s2  +  k  s  +  k  ) 
2  2  1 


and 


|G  (jto) 
2  1 


|G  (jto) 

2  3 


(1  +  Y  ) 

l 

2  ,  2  2~ 
(o  uo  +  y  ) 


2  2 
(Y  w)  +  (y^k  +  Y2) 


(19) 


(20) 


[to3  -  to(k  -  Y  k  )]2  +  [(k  -  Y  )“2  +  Y  k  j 

1  2  2  2  2  2  1 


But,  when  to=0,  ]G  (jto)  |  is  undefined  unless  Y  =  “l*  Under  this  condition, 
2 1  1 


we  now  have 


(to)  = 


2  2 

to  +  (y  -  k  ) 

2  2 


2  2  2  2 

[to  -  <o(k  -  Y  k  )]  +  {(k  -  Y)<0  +Y  k  ] 

1  2  2  2  2  2  1 


seT(<0) 


(21) 


Since  we  normalized  all  PDS's,  this  becomes 


(Y  k  )  s 

X'(to)  =  ■  fJ>  •  eT,(w) 


(22) 


(Y,  -  k  V 
2  2 

Let  a  =  [Y  ,  k  ,  k  ]  with  Y  =  -1,  then  compute  the  partial  derivatives  of 
,  —  2  12  1 

Se^t(to,a)  with  respect  to  a. 

DS  jCw,  a) 

__  _T _ 

All  the  information  is  now  ready  for  implementation  of  the  curve  fitting 

program.  The  criterion  function  to  be  minimized  is: 

CO  ,  .  2 

J(a)  =  /  1  [  S-  (co)  -  sg  t(to,  a)]  do) 

°  eT  T 

where  co  is  the  frequency  bandwidth  of  the  target  trajectory  (input)  and 

i  *■ 


(23) 


S- 


e^ (co)  is  the  normalized  PSD  of  the  sample  ensemble  mean  of  the  tracking  error. 


For  the  computer  program  we  used  the  discretized  form 

K 


'  o’  2 

J (a)  =  l  [  SeT(tok)  -Se  <cok,  a)]  Ato 
k=i  T 


(24) 


whare  Aco  is  the  chosen  increment  between  samples  in  the  frequency  domain  and 
0)^  =  k‘Aco,  co^  =  KAco.  The  iterative  equation  to  be  solved  by  the  program  is 


£ 


K 


3  .  ,  ,  ~  3  . 

-1+1  -1 


f  t 

.  -  p [  E  (  8  eT(V  -i)  )T  #  <  8SeT((\’  ~i)  ))"1 

k=l  «  *  a 

3  a  o  a 


K. 


<  E  (  aSeT(b)k’  V  )T  (  Se^(uk,  a1)  -  1  1 

k=l  3  a  r' 


I 

s- 


(25) 


where  p  is  a  step  size  factor  defined  in  Appendix  B.  The  values  of  model 

i 

parameters  obtained  through  the  identification  program  are  listed  in  Table  i 
below: 


AZIMUTH 

ELEVATION 

Y 

-1 

-1 

i 

Yz 

-1.77 

-1.877 

k 

25 

25 

1 

k 

2.6 

3.76 

2 

Table  1 

Once  these  gains  are  determined,  the  system  matrix  is  known.  The  second 
part  of  curve  fitting  involves  comparing  the  model  prediction  and  the  sample 
ensemble  standard  deviation  -5  D^,  of  the  tracking  error.  This  part  of  the  program 
was  computed  in  time  domain. 

To  begin,  we  need  the  solution  to  Eq.  (12).  (See  ref erence  [  11  ] . ) 

P(t)  =  <j>(t,  t  )  P(t  )  4>T(t,  t  )  +  <Kt, t)B  V(t)B  V(t,-r)dT  (26) 

ft  o  0  c  11 

0 

Select  P(t  )  to  be  a  diagonal  matrix  with  diagonal  elements  d  . ,  i  =  1,  2,  ‘ ' ' ,  n. 
o  -Li 

Thus,  the  second  diagonal  element  (t)  of  the  covariance  matrix  P(t)  is 
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(27) 


n  2  c  n  n  T 

P  (t)  =  Z  d..<p  .  (t-t  )  +  /  Z  Z  4>  .  (t-i)B  V(t)b  (t-t)dt 

22  i=l  11  2i  0  t  j=l i=l  21  1  1  2i 

0 

n  2  t  2 

=  Z  d  ,<p  .(t-t  )  +  /  (4>  (c-t)  +  <P  (t- t)]  V(T)dr 

i=^  il  21  0  to  22  2  4 

Since  K  and  F  are  known  and  since  G  (s)  and  G  (s)  can  be  obtained  from  Eq.(19), 

22  24 

we  need  only  take  the  inverse  Laplace  transform  of  these  quantities  to  obtain: 

<f>  (t)  =  eY2C 

22 


<P  (t)  -  -c  e  Y2C+c  e  C2tcos  c  t  +  c  e~C2tsi 


t  +  c  e  ^  sm  c  t 

3  4  3 


with  c.  defined  in  Table  2  below 

l 


Table  2 

V(t)  is  the  covariance  matrix  of  the  random  remnant.  It  has  been  found  that 
the  main  source  of  tracking  error  is  due  to  the  gunner's  uncertainty  about  the 
target  trajectory  dynamics  especially  the  target  velocity  0^  and  the  target  accel¬ 
eration  0j,.  Furthermore,  study  of  the  curve-fitting  between  the  empirical  data 
and  the  model  prediction  of  the  ensemble  standard  deviation  of  the  tracking  error 
has  indicated  that  the  gunner's  target  motion  uncertainty  is  the  dominant  part  of 
the  random  remnant  term  v(t).  Therefore,  it  is  proposed  that  the  covariance 
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function  V(t)  of  the  remnant  term  be  a  function  of  the  target  dynamics  as  follows: 


/*.  2  2 

V(t)  =  a  +  a  0  (t)  +  c  9  (t) 

1  2  T  3  T 


(29) 


where  a  ,  a  and  a  are  three  nonnegative  constants  to  be  identified  and  0T  and  8 
123  A  t 

are  estimated  target  angle  rate  and  acceleration,  respectively.  The  reason  that 
\  A  ^ 

-»  •  ** 

only  estimated  values  0T  and  8T  were  used  to  represent  target  dynamics  is  that  the 
gunner  doesn't  have  precise  information  about  ©T  and  0T  (i.e.,  target  uncertainty 
problem).  These  estimated  quantities  can  be  obtained  as  follows: 


\  =  U  )  -  (2  ) 

1  “11  ~13 


(30) 


where  z_  is  the  state  of  the  closed  loop  AAA  system  with  the  observer  model,  and 
(z  )  and  (z  )  are  the  first  and  third  components  of  the  vecor  z  .  Then  by 

—li  —13  —l 

approximation: 


%(tk)  =  tT(tk)  -  e^Vi) 


At 

where  At  is  the  sampling  interval,  t^  =  k  •  At.  Let 

b  =  [a  a  a  ] ,  then  in  addition  to  P  (t,  b),  we  need 


l  2  3 


22 


9  P  (t,b) 

22  — 

8b 


(31) 


Finally,  the  criterion  function: 

l  t 


f 


J  (b)  -  /  [  SDft)  -  P  (t,  b)  ]  dt 

—  „  T  22  - 


(32) 


For  the  computer  program,  a  discretized  form  of  the  criterion  function  is  used. 


J  (b)  =  I  [  SD‘  (t.  )  -  P  (t,  ,  b)  ]2  •  At 


k=l 


T  '  k 


22  k  ~ 


(33) 


The  iterative  equation  to  be  solved  by  the  program  is: 


b...  =  b.-p 
— i+1  -i 


r  /9P  (t  ,b.)  T 
L  (  22  k  —L  ) 

9P  (t  ,b.) 

(  22  k  — 1  ) 

-k=1  % 

8b 

K 

£ 

k=l 


E  (3P22(tk’^L))T  (P  (t  ,b  )-SD  (t  ) 

-  •  22  k  —i  Ik. 

9b 


The  values  of  model  parameters  obtained  through  the  identification  program 
are  listed  in  Table  3. 
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AZIMUTH 

ELEVATION 

a 

.001 

.0005 

1 

a 

.025 

.05 

2 

a 

.008 

.0025 

3 

Table  3 


A  block  diagram  and  computer  program  listings  for  each  of  the  two  parameter 
identification  parts  (i.e.,  ensemble  mean  and  ensemble  standard  deviation)  can 
be  found  in  Appendix  C. 
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SECTION  IV 


■  ‘  '  —  ‘  __  -|f|  ,  |  -  .  -||  |  1  ;  r-  .  — ... 


COMPUTER  SIMULATION  RESULTS 

A.  Discretization  of  Observer  Model  Equations 

To  more  efficiently  use  computer  time  and  core,  the  observer  model 
equations  were  first  discretized  and  then  programmed. 

So,  given  the  continuous  state  variable  equation  (10),  discretization 
yields: 


^n+1 


A  At  ,  rAt  Ai(At-s)„  o  .At  A! (At-s)_  . 

=  e  1  z  +  /  e  F  0„(t  +  s)ds  +  /  e  ■  B  V(t  +s)ds  (34) 

o  l  T  n  o  in 


— n 


=  eAlAtz  +  /  C  eAl(At_s)ds  F  6  (t  )  +  /  eAl  (At“s)ds  B  v(t  ) 
— n  o  i  T  n  o  in 


=  eAlAtz  +  /  eAlU  da  F  0  (t  )  +  /  e^10  da  B  v(t  ) 


At 


AjO 


At 


Aia 


— n 


i  Tv  n  o 

Hence,  we  have  the  following  difference  equation: 


l  n 


z  . ,  =  $  z  +  T  6  ,  +  T  v 

— n+1  —a  i  1  n  2  n 


(35) 


where 


.  AiAt  r  .At  AiO,  .  v  p  _  ,AC  A.o 
$  =  e  ,P  -  J  e  da  F  *  T  /  e  1  da  B 
10  12  0  1 

At  is  the  sampling  period,  0_  =  0_(t  ),  and  v  is  a  random  sequence  with 

1  ,n  in  n 

the  following  properties: 

E  (v  ]  =  0 


E  [  (v  -E[v  ])(v  -E[v  ])T]  =  V  = 

n  n  n  n  n  At 


(36) 


where  V(t)  is  defined  in  Eq.  (7). 

Taking  expectation  value  of  both  sides  of  Eq.  (35).  We  get 

z  =$z  +  r  8 

-n+1  -n  i  T,n 


where 


Ai+i  *  EIVi!' 


(37) 


The  covariance  of  z  is  defined  as: 

—n+1 

X„M  “  EI(Vl  -  W  (Vl  -  Vl>T'- 

It  can  be  shown  [11]  that  is  governed  by  the  following  equation 


_ _ f 


20 


x  ,  ,  =  $  x  <tT+rvrT 

n+l 


(38) 


n  2  n  2 

The  predicted  ensemble  mean  tracking  error  can  bs  obtained  from  the  second 
element  of  of  Eq.  (37),  and  the  predicted  covariance  of  the  tracking 

error  can  be  found  in  the  second  diagonal  element  of  the  covariance  matrix 
Xn+1  of  Eq.  (38). 

B.  Simulation  Results 

The  numerical  values  of  the  parameters  of  this  gunner  model  were  determined 
in  section  III  with  respect  to  the  gunsight  dynamics  system  (Eq.  (1))  and  a 
deterministic'  target  trajectory.  Angular  information  from  Trajectory  4  was 
the  input  to  the  curve-fitting  program.  Since  the  parameters  have  been  selected, 
all  the  necessary  matrices  in  Section  IV  A.  are  defined.  So  the  AAA  gunner 
model  is  now  ready  to  be  used  for  computer  sumulation.  A  computer  program,  OMS, 
for  simulating  an  AAA  system  with  this  model  representing  the  gunner  response 
was  implemented.  A  block  diagram  and  program  listing  of  this  simulation  can  be 
found  in  Appendix  D.  The  input  to  this  program  is  the  trajectory  of  the  target 
motion.  Initially,  only  Trajectory  4  was  used.  The  outputs  are  the  model 
predictions  of  the  ensemble  mean  and  standard  deviation  of  the  tracking  errors. 
The  results  for  Trajectory  4  are  plotted  in  Figures  6a  through  6d.  Each  graph 
contains  empirical  data,  observer  model  predictions  and  P001  formula.  (P001 
formula  refers  to  a  simple  formula  to  predict  tracking  error  used  in  the  P001 
attrition  model  program  [12].) 

Azimuth  mean  tracking  error  and  standard  deviation  are  shown  in  Figures  6a 
and  6b  respectively.  Similarly  Figures  6c  and  6d  show  results  for  elevation 
tracking  errors.  It  is  obvious  that  the  predictions  by  the  P001  ormula  did 
not  match  empirical  data  well.  However,  matching  between  empirical  data  and 
observer  model  predictions  was  very  good.  All  these  results  indicated  that 
this  AAA  gunner  model  is  able  to  represent  the  trend  of  the  gunner  response  in 
the  tracking  task.  It  was  also  noted  that  the  sharp  peaks  in  the  empirical 
data  curves  were  not  matched  by  the  model  predictions.  This  may  be  due  to 
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the  simplified  gunsight  dynamics  used  in  designing  this  model  or  an  insufficient 
number  of  runs  used  in  generating  the  sample  ensemble  mean  data.  Further  re¬ 
search  concerning  these  problems  is  necessary.  Next,  this  gunner  model,  with 
the  same  set  of  parameter  values,  was  used  in  OMS  to  predict  the  tracking 
errors  for  three  other  target  trajectories  (1,2,3,  in  Figure  2).  Figures  7,  8, 

9  picture  these  results  in  the  same  format  as  listed  above.  All  the  simulation 
results  show  that  this  AAA  gunner  model  with  the  same  parameter  values  gives 
model  predictions  in  good  agreement  with  empirical  data.  Therefore,  the  observer 
model  is  a  predictive  model  in  the  sense  that  it  can  be  used  to  predict  tracking 
errors  of  an  AAA  system  for  various  target  trajectories  with  similar  frequency 
bandwidths.  In  addition,  the  observer  model  is  also  an  adaptive  model  since  its 
parameters  depend  on  the  gunsight  dynamics  and  the  target  trajectory. 

A  comparison  of  the  model  prediction  accuracy  between  the  observer  mode.' 
and  the  optimal  control  model  has  been  done  for  these  4  target  trajectories. 

All  the  results  showed  that  both  models  give  accurate  predictions  of  the  tracking 
errors.  Figures  10a  through  lOd  show  the  ensemble  mean  and  standard  deviation 
of  the  tracking  error  as  predicted  by  the  optimal  control  model  for  both  azimuth 
and  elevation  tracking  task.  Upon  comparing  Figures  6  and  10,  it  is  obvious 
that  the  AAA  gunner  model  (i.e.,  observer  model)  developed  in  this  paper  can 
predict  the  tracking  errors  as  accurately  as  those  by  the  optimal  control  model. 
The  computer  simulation  time  of  the  closed  loop  AAA  system  using  the  observer 
model  is  only  6.5  seconds,  while  37  seconds  of  simulation  time  are  needed  to 
execute  the  optimal  control  model.  Therefore,  a  reduction  of  85%  computer  simu¬ 
lation  time  can  be  obtained  by  using  the  observer  model  instead  of  the  optimal 
control  model.  So  this  AAA  gunner  model  based  on  the  observer  theory  is  very 
useful  in  the  analysis  of  the  performance  of  the  AAA  gun  system. 


SECTION  V 


CONCLUSION 

The  Luenberger  observer  theory  has  been  applied  to  design  an  AAA  gunner 
model.  The  key  design  requirement  was  to  make  the  model  structure  simple  so 
that  it  needs  much  less  computer  simulation  time  than  the  other  models.  It  was 
also  important  to  predict  the  tracking  error  accurately.  Both  specifications 
have  been  met.  In  addition,  this  report  has  presented  a  parameter  identification 
program  which  can  easily  determine  the  numerical  values  of  the  parameters  of 
this  model.  Then  the  model  is  ready  for  the  computer  simulation  of  the  AAA  gun 
system.  The  Aerospace  Medical  Research  Laboratory,  WPAFB,  has  applied  this  model 
to  the  study  of  a  foreign  AAA  gun  system.  This  model  is  now  being  applied  to 
study  several  other  foreign  AAA  gun  systems  and  SAMs.  The  identity  observer 
used  in  the  observer  model  still  possesses  a  certain  degree  of  redundancy.  The 
redundancy  stems  from  the  fact  that  the  identity  observer  approximately  con¬ 
structs  an  estimate  of  the  entire  state  but  part  of  the  state,  as  given  by  the 
system  outputs,  are  already  available  for  direct  measurement.  This  redundancy 
can  be  eliminated  by  the  use  of  a  reduced-order  observer  to  replace  the  identity 
observer.  Further  research  is  continuing  to  develop  an  observer  model  which  uses 
fewer  states.  This  will  further  simplify  the  structure  of  the  current  observer 
model  and  shorten  the  computer  time.  Now,  to  determine  the  parameter  values 
for  the  model,  it  is  necessary  to  have  empirical  tracking  data  from  a  given 
weapon  system  available.  So  another  worthwhile  extension  of  this  work  is  to 
develop  parameter  adjusting  rules  for  this  model  such  that  the  numerical  values 
of  the  parameters  can  be  obtained  without  the  use  of  the  empirical  data.  This 
project  is  now  under  study. 
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FIGURE  10B:  OPTIMAL  CONTROL  MODEL  PREDICTION  (PHI  FLYBY) 


FIGURE  10D:  OPTIMAL  CONTROL  MODEL  PREDICTION 


APPENDIX  A 


A  Comparison  of  the  Human  Operator  Models  for  AAA  Systems 

This  appendix  will  give  a  general  description  and  comparison  to  the 
following  models: 

1.  McRuer  Crossover  Model  (based  on  classical  control  theory) 

2.  Optimal  Control  Model  (based  on  optimal  control  and  estimation  theory) 

3.  PID  Structure  Modified  Optimal  Control  Model  (simplified  Optimal  Control 
Model) 

4.  Observer  Model  (based  on  the  observer  theory) 

1.  McRuer  Crossover  Model 

The  structure  of  the  McRuer  crossover  model  is  described  in  the  frequency 

domain  in  Figure  11(a).  yq(s)  is  the  transfer  function  of  the  controlled  system 

and  the  crossover  model  is  composed  of  a  linear  element  Y  (s)  and  an  random  element 

n 

(remnant).  An  interesting  relationship  of  the  open  loop  transfer  function  Y  (s)*Y  (s) 
was  found  by  McRuer  as  follows: 

Yh  (ju>)  Yc(j<0)  S  v"jt0Te 

jw 

around  the  region  of  crossover  frequency.  The  two  parameters co^  and  are  crossover 
frequency  and  effective  time  delay  respectively. 

This  model  was  developed  using  classical  control  theory,  i.e. ,  frequency 
domain  analysis.  It  is  a  useful  model  of  a  human  operator  primarily  applied  to 
single-input  single-output  linear  time- invariant  systems  with  random  or  random¬ 
appearing  inputs.  The  applications  to  multivariable  systems  or  time-varying 
systems  or  deterministic  input  forcing  functions  are  not  straightforward. 

2.  Optimal  Control  Model 

Optimal  control  model  assumes  that  the  human  operator  behaves  in  an  optimal 
manner  subject  to  his  inherent  limitations  and  to  the  requirements  of  the  control 
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task.  Kalman  filter  is  used  in  the  model  to  obtain  the  optimal  estimation  of  the 
state  of  the  dynair  _ystem.  Then  the  optimal,  linear,  state  variable  feedback 
control  is  implemented  by  using  the  estimated  state.  The  block  diagram  of  the 
model  is  shown  in  Figure  11(b).  Optimal  control  model  is  b^sed  on  the  state 
space  techniques  to  describe  the  dynamic  system  and  hence  is  convenient  to  be 
applied  to  multivariable  systems  or  time-varying  systems.  So  the  optimal  control 
model  is  a  good  gunner  model  for  complex  weapon  systems.  The  problems  of  the 
optimal  control  model  are  summarized  as  follows: 

a)  The  structure  of  the  optima]  control  model  is  complicated.  There  are 
two  Riccati  equations  (which  are  nonlinear  matrix  differential  equations) 
involved  in  the  model  to  be  solved.  Therefore  the  computer  execution 
time  to  generate  model  predictions  of  the  tracking  errors  using  optimal 
control  model  is  lengthy. 

b)  There  is  no  systematic  method  to  determine  the  parameters  of  the  optimal 
control  model.  Usually  it  takes  a  long  time  to  modify  the  cost  function 
or  model  structure  in  order  to  obtain  a  set  of  parameter  values  to  give 
satisfactory  model  predictions.  Generally  speaking,  these  parameters 
are  determined  through  trial  and  error. 

3.  PID  Structure  Modified  Optimal  Control  Model 

The  PID  structure  modified  optimal  control  model  is  a  simplified  version  of 
the  optimal  control  model.  Its  structure  is  shown  in  Figure  11(c).  Obviously,  it 
is  much  simpler  than  the  optimal  control  model.  In  addition,  this  model  assumes 
thatthe  human  operator  learns  an  input-output  type  of  internal  model  relating  the 
human's  control  u(t)  to  the  displayed  variable  e^(t)  (i.e.,  input  to  the  human 
operator),  without  using  any  explicit  representation  of  the  actual  system  dynamics 
and  forcing  functior.  The  advantage  of  this  idea  is  to  have  simple  mathematical 
equations  for  the  model.  The  disadvantage  of  this  formulation  is  that  it  may  be 
difficult  to  tune  the  parameters  effectively  to  obtain  accurate  model  predictions. 

Another  good  idea  in  this  model  is  the  consideration  of  the  lower  order 
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internal  model  of  the  controlled  plant.  But  this  causes  the  divergence  of  the 
covariance  matrix  of  the  Kalman  filter.  Further  study  is  necessary  in  order  to 
use  this  idea  successfully. 

4.  Observer  Model 

Two  of  the  main  design  requirements  in  the  development  of  the  observer 
model  are: 

a)  systematic  determination  of  parameters  of  the  human  operator  model, 

b)  simple  structure  of  the  model  (this  shortens  execution  time  in  the 
computer  simulation  of  an  AAA  tracking  task). 

The  Luenberger  observer  theory  was  applied  to  develop  this  model  and  its 
basic  structure  is  shown  in  Figure  11(d)  for  comparison.  In  this  figure,  the 
observer  plays  the  role  of  the  Kalman  filter  and  the  feedback  controller  plays  the 
role  of  the  optimal  controller.  An  observer  generates  an  estimate  of  the  state 
for  a  dynamic  system.  Its  advantage  is  that  no  Riccati  equation  is  required  to 
compute  the  observer  gain.  (pote  that  in  the  optimal  control  model  to  compute  the 
optimal  Kalman  gain,  a  nonlinear  matrix  differential  equation  has  to  be  solved.) 
Similarly  there  is  no  quadratic  cost  function  involved  in  the  formulation  of 
observer  model.  Only  a  linear  state-variable  feedback  control  law  is  used  in 
this  model.  Hence,  the  structure  of  the  observer  model  is  much  simpler  than  that 
of  the  optimal  control  model. 

In  addition,  the  parameters  of  the  observer  model  are  not  determined  by 
trial  and  error  tuning.  The  developed  parameter  identification  program  (i.e., 
least  squares  curve-fitting  program)  in  Section  III  can  identify  the  parameters 
systematically. 
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APPENDIX  B 


Development  of  Gauss-Newton  Gradient  Method 


Consider  a  scalar  function  J(a)  of  a  vector  parameter  a  and  its  first 
derivative  3J(a)/3a.  Taking  first  order  Taylor  series  expansion  of  3j(a)/3a  with 


respect  to  a  minimum  a*  of  Jfe),  we  have: 

2 

3J (a)  3J (a*)  3  j(a*) 


[a  -  a*] 


(Bl) 


3  a  3  a_  3a3^ 

The  first  term  of  the  right  hand  side  of  Eq.  (Bl)  is  zero  because  a*  is  a  minimum 
of  J(a).  Now,  if  the  second  derivative  of  J(a)  with  respect  to  a  at  a*  is  not 


singular,  then  Eq.  (Bl)  can  be  rewritten  as  follows: 
2 


3  J(a*) 

1*  -  a  -  [ - =H  .  ( 


3JCa) 


) 


(B2) 


3a3ci  3a_ 

Eq.  (B2)  can  be  used  as  a  basis  for  an  iterative  relation  to  obtain  the  minimum 


of  the  function  J(a).  The  iterative  relation  is  given  by 

2  -1 


a  =  a 

-i-rl  % 


3  J(a..) 
3a3aT 


3J(a.) 

3a 


(B3) 


This  is  the  Newton  gradiant  method.  In  order  to  satisfy  the  descent  property, 


i.e.,  J (ai+^)  <  J (a^) ,  a  one-dimensional  search  for  is  performed, to  adjust  the 
step  size  of  increment  in  the  iterative  procedure.  Newton  method  can  now  be 
written  as:  „  -1 


a.  -  p. 
1  x 


'  32J(a.)‘ 


3a3a 


3J(a.) 


3a 


(B4) 


The  Gauss-Newton  technique  is  an  attempt  to  obtain  convergence  characteristics 
that  are  similar  to  the  Newton  method  without  calculating  the  second  derivative 
of  the  criterion  function  with  respect  to  a.  Consider  a  known  function  f^(x)  and 
a  second  function  f(x,£)  where  a_  is  a  parameter  vector  of  dimension  p,  to  be 
determined  in  order  to  minimize  the  criterion  function: 
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J(a)  =  /  [f1(x)  -  f(x,  a)  ]  dx  (B5) 

xo 

Where  x0,  x  are  integral  limits.  Take  a  first  order  Taylor  series  expansion  of 
f(x,  a)  about  a  certain  initial  guess  £0  and  substitute  into  J (a) . 

J(a)  =  /  [  f  (x)  -  f(x,  a0)  -  -o^  (a  -  aQ)  ]  dx  (B6) 


For  simplicity  we  will  assume  that  f^(x)  and  f(x,  a)  are  acalar  functions  of  x 
and  (x,  a)  respectively,  and  a  is  a  p  dimensional  vector.  Now  the  scalar  criterion 


function  J(a)  can  be  rewritten  as: 


J  (a)  =  /  { (f  1  (x)-f  (x  ,a  ))2-2(f  (x)-f(x,a  ))  (a  -  a.)  +  - o •* 

—  i  —  o  1  —  o  -  —  — u  -  va  -  «</ 

o  3a  [_  3a 


Taking  the  partial  derivative  of  J  (a)  with  respect  to  a_,  we  get 


f  -2(f1(x)  -  f(x,  a0)  )  - 

xo  8a 


3f(x,  a0)  3f(x,  a0)  3f(x,  a^) 

- +  2  -  (a-a0) -  dx 


Note  that  the  gradient  of  the  criterion  function  is  zero  at  its  minimum  a*. 


3J (a*)  Xf 


3f(x,  a0)  3f(x,  a  ) 


/  ~2  (f1(x)-f(x,  a0)) 

xo 


a  )  3f(x,a0) 

-  (a*-a  ) - dx 


3a  J 


So  taking  the  transposition  of  the  resulting  equation  gives  us 

X£  3f(x,  a  )  T  Xf  3f(x,  a0)  T  3f(x,  afl) 

f  (  - )  (f,(x)  -  f(x,  aQ)  )dx  =  /  ( - )  -  (a*-a0)  dx 

xn  3a  xn  3a  3a 


For  the  discrete  time  case,  the  analogous  result  is: 


K  9f(xv.a0^T 

Z  ( - ± - ) 

k=l  3a 


T  K  3f(x  ,  afl)  T  3f(x  aQ) 

.)  (f  (x  )  -  f(:c  ,  a0)  )  =  I  (  - - -  )  (  - - - )(a*-afl) 

k=l  3a  3a 
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Hence , 


-1 


a*  =  *0  + 

K  /3E(xk„0)V  3f(xk,a0) 

fl(xk)-f(xk,  a0) 

k=l\  3a  3a_ 

k_1  \  da 

The  iterative  relation  is  given  by: 


(B13) 


In  order  to  satisfy  the  descent  property,  i.e.,  JCa^^)  <  a  one  dimensional 

t 

sear-.h  for  p  is  required.  Thus  the  modified  Gauss-Newton  iterative  relationship 
is  expressed  as:  ~1 


a  1  =  a.  -  p 
—l+l  —i 


-K  /af(xk,ai)ypf(xk,a1)V 
k=l  \  3a  '  '  3a  ' 


^  /3f(xk»^AT 


k=l 


da 


)  (f(xk’V-fi 


(xk) 


(B14) 


Let  the  second  term  of  the  right  hand  side  of  Eq.  (B14)  be  denoted  by  +  pD. 

— i 

then: 


a  =  a,  +  PD. 
—l+l  -i  -a 


Initially  p  =  1.  If  the  descent  property  is  not  satisfied  chan  p  is  halved, 
process  continues  until  J(a_+^)  <  J(a^)  or  a  iow^r  limit  is  reached.  At  each 
iteration  p  is  reset  and  the  criterion  function  is  again  checked. 

Let 


(D.). 

E-  = 

3  (a). 

-i  1 

where  (D . ) .  and  (a, 

-l  j  -1 

convergence  of  Eq. 

lEil  < e 


for  j  =  1,2,. . . .  p. 
til 

)  are  the  j  components  of  and  a^  respectively. 
(B14)  is  established  whenever 
for  j  =  1,2, - p. 


The 


This 


where  G  is  a  preassigned  small  positive  number.  Otherwise,  further  iteration  is 
required. 

Following  is  a  flowchart  for  a  computer  implementation  of  the  modified 


Gauss-Newton  method. 


(B12) 


) 
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APPENDIX  C  -  PARAMETER  IDENTIFICATION  PROGRAM 

A.  FLOW  CHART  OF  PARAMETER  IDENTIFICATION  FOR  PARAMETERS  ASSOCIATED  WITH  MFAN 
TRACKING  ERROR  EQUATION 


READ : 

F 1  (I)  =S-  (o') 

eT  1=1,..., N 
F2(I)=S„  (to) 

T 

aQ(k)=  rr2,k1,k2  ] 
k=l,2,3 


INITIALIZATION?  " 

F3o=s^T2y(w»^> 

Bo=  ^SfjXZT^’^ o) 

~  ^aOc) 

N 

J0(-)=A(F1(I)-F30(I)): 


CALCULATE 


Fi+l=Vp£l 


UPDATE:  F3i+1 


■Hi/  Ji+1 <  Ji 


P  =  p/?  L.T^S/^  CHECK:  \  IQ  /•  stop 

\  p  <  limi  /  Terror  i 


COMMENTS:  EXECUTION  OF  PROGRAM  FIT1 


1.  COMPILE  PROGRAM 

2.  INPUT:  ATTACH  TAPE1  AND  TAPE2  where 

TAPE1  CONTAINS  EMPIRICAL  STANDARD  DEVIATION  OF  TRACKING  ERROR 
TAPE2  HAS  PARTIAL  DERIVATIVE  INFORMATION  (described  later  in  SDCON) 

3.  ATTACH  D.L.  KLEINMAN  [13]  LIBRARY 

SUBROUTINE  -  GMINV ,  VMATI,  VSCALE 

4.  RUN 

5.  OUTPUT:  PARAMETERS  VALUES  FOR  [dp  ct2>  a^] 


X?  MPUT  >  OUTPUT, TA*£  l , TAP£2  > 

-  -iKS-.Sj  •  1  ~°  G0  TV°^  0  TO  STOP 
..  .  1f<TFG,EQ,0)STO° 

Ream  i 
REWIND  2 
CALL  INI7 
GO  TO  l 
tNO 

subroutine  InIT 
INITIALIZATION  ROUTINE 
k£AD  IN  OATA  TO  3:  PTTT-n 
MAKE  INITIAL  r,u-sI 
COM*ON/VAC/LTmi,L  H 

C0HH0N/NAT/FHH2  =£7 I , 

OATA  N,H,  13/112-:,  .'34  i/  °  ’  3ai2E)tA(31,AA(3>fB(U25,3},0(3> 

PPINT* *  30  H  Tror  j'l  LI^t.LlN’ 

REAO*«LI.N  1,LIM2,I"  ’  1 

00  10  1*1,1 

READ (1,31 Z,Zi,Fl  (T)  ,7? 

Fl(ItsPi(T)*Cl(I) 

FORMAT  (<*Gl2.5) 

N0£*:  “"“f*  8  18  C0BPut*d  *»  *  separate 

CONTTm.  r  r0Utlne  SDC0N  (Hating  folios) 

PRINT*  r-  r  -  f^nCe  B  refflaljls  constant  throughout 

REAO*,AA  '  *  *  i  »UtSS— Aw3  2.  A  ^h£,  <f£SPtion  Process. 

PPINT*,  1  OH  nr  TU-S3,AA 
0£L*l 

CALL  COE*  <U,I°) 

CALL  JNlNnjJj 

PRINT*, 10k  (ej? 

CALL  LOOP  (*J2, 1=) 
tND 

subroutine  jn:n(ej) 

SuSI,/MftVWlfU2n,P?ai2,',,r5ai25,,fl‘3».AA(3),  3CU25, 3), 0(3, 

00  13  I=i,i 

SUM’SUN*  <  FI  (11  -ft  (!,,**■> 

RJ=SUM 
END 

SUBROUTINE  00(1=) 

CALCULATE  3  NATRIX  «CN  0  AND  » 

0IM-.V5I0N  R<31 ,0(  3,  3,  ,H1  (3,  ?; 

COMMON/VAR/LIMi,Li;?;~™;fTL 

CONMON/NAm/NOIM‘,NOT^lic!!r3(wS’.fl<;!,’AAU,’e<1125*5»  =0  C3I 
NCIM=IS3Nr  <l-T=*  , ,",1-',1»>0«*9t5)/lNeO/<IN,KO«T 
KINsKOUT**' 

00  13  I*1,IP 

failiS ^  .  •.  i.  -  iww* iw*** 

00  35  «■!.«  *£  JS  **%  « »e 
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SN=F3 (K) -  FI (<) 

CC  30  1*1, IP 
£.0  :>  j=i,:o 

CCI,  J»:0€],JWO(K,I|*!i(K,  ji 
«<I)  =  S.<I)  ♦1K,!)*S'I 
CONTINUE 

CALL  G'HNVt:3,I»,0,wlt^t0) 

CALL  VMATl(Ml,!»,i«;,  io.oj 

call  vscAL-:n,o,r=,-i,> 

S.NO 

SU3R3UTINE  wD3«fOj?,XP| 

ITIRA^IDM  >,0C;3S 

COM  PUT  C  A  C*1 )  =1  ( I)  O_L*0  (I ) 

OJTHMINC  «HT*.  A (1*1)  IS  A'ok  pTjpi 
DIMENSION  100 (-1 
COMMON/tfAS/LIMl,LM?,C--,H.  l.OiL 

N??^{=5.1,CUU?e,*r?(ll25,’r3<1125>’A,3>’AM3,’8(l^^>‘0<3» 
CALL  03 ( I rl 
00  1 3  C  1=1,13 

om*C(D  «)-:l 

m (I ) =AA( I  )♦!(!) 

CAL v  COIF  n» 

CALL  JMINnjl) 

1F(PJ1,LT , ?J2>SC  To  30 

ofl*o:l/2, 

rtCT  sMC'*  1 

iF(NCT.Ll,LlHi)SO  TO  ? 

?mr..rMS«3P  1.3H  ,CT-.,OT.;h  .MCT=,M;t,pji,RJ2 

oo  3;  :=i  ,:3 
b00(I)s)3£(3rr)/AA(I)) 

I F ( ooi ( i > , ST.- "160  *o  Pj 

GC  TO  -o 

MCT=irT*l 

.lr<MC7,Li.LlM:)S3  To  c 

PRINT*, 7H£, 2,Jm  4=,fl,:s  tl-Ts,,CT,-H  hCT.,NCT,PJ1,RJ2 

RJ<!  =  ’J1 
NCT  =  0 
DEL=1, 

CO  25  1  =  1, 1» 
mA(IJ=6(I) 

GO  TO  l 

PRINT*, 4  =  ,A.-‘-‘|-Ts,»4OT,.HMCr=.*0T 

.'’«INT*,3HM=.A4 

PRINT*, 3H  3s, o 

PclNT^,RHwlT'l*,'Jl 

lN3 

5U3R3UT:n=  03"e«H,IP) 

CCMdJTi  FUOTION  OIVLN  ’A- AMt  T  £RS  A 

PAP.TIflL  0-,iWA"lv.S  vr  -7  ( 4 ,  A  )  WCT  A  nR-.  In  MAT =  I X  0 
OlMiMSIDN  A(IP)  A  ' 

COHNON/tfi  b,LI'11.LT',,,7,:.h.n,C_L 

GO^IO  1=1  |^ri  (11?-  *  •"*.(  11?^»  ,r3  (112;1 ,  A  (3) ,  AA(3>  ,5  (1125,3)  ,0(3) 

CQNT'NUT  *3**<*^’***W^-*+3(I,2)*N(?>*2(it?)'*N|3) 
l  NO  * 


„  -am-*— 

raoa 


B.  FLOW  CHART  OF  PARAMETER  IDENTIFICATION  FOR  PARAMETERS  ASSOCIATED  WITH  STANDARD 
DEVIATION  OF  TRACKING  ERROR  EQUATION 
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COMMENTS:  EXECUTION  OF  PROGRAM  FIT 

1.  COMPILE  PROGRAM 

2.  INPUT:  ATTACK,  TAPE1  and  TAPE 2  where 

TAPE1  CONTAINS  PSD  OF  EMPIRICAL  MEAN  TRACKING  ERROR 

TAPE2  CONTAINS  PSD  OF  6  FOR  GIVEN  TRAJECTORY 

T 

3.  ATTACH  D.L.  KLEINMAN  £1 3 ]  LIBRARY 

SUBROUTINE  USED  -  GMINV ,  VMATl,  VSCALE 

4.  RUN 

5.  OUTPUT:  PARAMETER  VALUES  FOR  [Y  2,  kX>  k2  ^ 


PROG^M  F  it  <  I  N'°U  TVO  U  TP'uf ,  TAP  £1 ,  T  A  Pi  2 » T  AP  E  6 ) 

PRINT*, 50HTYPI  1  TO  GO  OR  TYP£  0  TO  STOP 

■REAO^IF'G -  “ 

IF(IFG,£Qj3>ST03  _  _  _ 

REWIND  l 
REWIND  2 
CALL  iNIT 
GO  TO  1 

END . . 

SUBROUTINE  INIT 

INITIALIZATION  ROUTINE 

READ  IN  DATA  TO  31  cITTEO 

MAKE  INITIAL  GUESS 

C0MN0N/VAR/LtMl«LIM2,cE,IP,H.»N,0!:L 

C0M0N/NAT^Fr(5i3f,P2(5T3l,F2(5l3)  ,A(3),AA(3)  ,3  (£'13, 3) ’,D  (TF  ”* 
DATA  N,H,  1= /51 3* «  G2-,**l»3/ 

PRINT*, 30H  TY°£  IN  LIM1 , L IH2, EE  ’ 

REA O*, LIN  1,  LINE, \  • 

PRINT  *,3H  4=,N,2HH=tW,5HLlU=tLIKl,5HLlN2=,LIM2,3Ht£=tc*,3HlP 
00  10  1=1, N  ___ 

READ(1,3) Z,Pl(I> .  ’ 

RLA3(2,h) ZtF2(I) 

FORMAT  (IK  ,2012,-,) 

FORMAT (IX  ,2G12»5> 

CONTINUE 

PRINT*, 30h  FIRST  GUfSS 

RtAO* , A  A  . ~ .  . 

PRINT* ,  1 0  P  137  GUESS, AA 
DEL=1 

CALL  COIF  (AA) 

CALL  JNIN  (RJ2) 

CALL  LOO»(<J2) 

FPE  0=5,  " . ~ .  * . 

DO  20  1=1,1 

WRITE  (6,2  5)  FRIO,  F3  (I) 

FP.EQ=FRE0+1 

FORMAT (2G 12, 3) 

END 

SUBROUTINE  JNIN(RJ)  .  ‘  . 

compute  5 j=crit:rton  function 
C0MM0N/VAR/cINl,LIN2,EE,IF,H,N,D'L 

C0NM0N/^AT/F1(513),F2(513) ,F3 (513) , A (3) ,AA<3> ,3(513,3) ,0(3) 
3UM=3. 

CO  10  1=1 ,1 

sum=3um+  (Ft(i)-F*i  mj*»r " "  "  . . 

RJsS'JN 

END 

SU3R0Ui'In:  00 

CALCULATE  0  MATRIX  FROM  0  mNO  R 
CIMENSION  R(3),0(3,3),W1(?,3) 
COMMON/VAA/L'!'fl,LIN2,:£,IP,N,N;OEL 

COMMON/MA VF1 (513) ,F?(=13),F3 (313) ,A(3) , AA(3) ,3(513,3) ,0(3) 
COMMON/MA  I (1/NOIH ,Nr)IPl',COMi  ( 3,  T)  /INOU/KIN,  KOUT 
N0IM=3  jND  INI8** 

KIN=KC'JT  =  6  ^ 

00  10  1=1,1=  TV*** 

R(H  =  0.  . . 
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Q  (I , J)  =  0. 
uC  li  <*1,4 
3N=r3 (<) *  C1 (<) 

CO  13  1  =  1 
DO  25  J= 1  .1 3 
2^  G(I«J)=Q<  I.JJ»3(K,I>*'5C'<,  J) 

30  %(I)=K(I)  ♦1«.I)*SN 

35  wCNTINJ^ 

CALL  G'HNV<I=M3,n,Wi,M^,0» 
cALL  V9A71(W1,*.IP,IP,0) 

CALL  V3CAi:tD»3.ID.-l*f 
c.N0 

iUlROUTIU;  *.30®(5J2) 

C  ITERATION  7R00i3S 

C  CCMPJT-.  aCHI^Ct+DiL^CI 

C  Cs.Tc.RMIN*  4H2N  *  C 1+ 1  >  13  A  3Ci  PTAGci 

i :m£nsion  3oon» 

COM'tON/;Afi/Lm.LH2.r*,  IP,H,n,02L 

cOMMON/MA  7/Fl  (51?)»c2(3l3)»r3(5l?)»A(3).AA(3)  ,3(513,3), 0(1) 

NCT  =MC r=3  < 

1  CALL  DO 

2  JO  130  1*1.1* 

3(i)oo  *i:l 

100  A<I>=AA<I)0(I) 

CALL  CO=P(\) 

5  CALL  JMINRJ1) 

6  1F(RJ1.L7.5.J2)G3  T0  33 
J-L=3LL/?. 

NC7  =NC7* 1 

IF(MCT.L-.lIM1)G3  TO  2 

PRINT*, 7M£R->DR  1*  3H  A=,4t>H  N0T=,N07 ,5H  MCT*,HCT,R J1,*J2 
STOP 

30  00  3:  1  =  1  ,t3  . 

000(I)=ABSO(I)/1ACI) ) 

35  IF(D30(I)  .Sr,;i)GO  T(?-'23  ”  "" 

GO  TO  *0 

20  MCT=1CT»1  . ~  ' 

if(m:t.l£.l:i2)go  to  s 

PRINT*,7H£RRQ?  2,  1H  A=  ,‘A f  -I NC'f  =,NCf  , pH  MCT= ,MOT , *J1 ,RJ2 
STOP 

9  RJ2=RJ1 . 

NCT  =  0  _ 

0EL=1.  . 

00  25  1=1  ,1° 

25  AA ( I) =A (I ) 

GO  TO  1 

40  PAf NT »o, TA'(I)Ti=l ,  f6),NCf  i  iZf  . . . . 

PRINT*, 3HAA=,AA 

PFIMT*,3H  3=V3  . "" 

PRI  NT  * ,  5  H  JHN=,RJ1 _ F>tt,  Ts  BIST  *W,m  WJgW" 

45  FCRH4f(iX,3G12.5.2X,2I5)  '  Pt.  i  ^JtiSD  WSM  ^  * 
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tN0  _ 

'  ■■SUBROUTfN-  COrr  (M>  . .  . 

C  COMPUTE  FACTION  F3  (W,A_)  GIVEN  PARAMETERS  A 

C  EXPRESS  PARTIAL  OIRiVATIVt S  Oc  F?IW,A>  WRT  A  IN  MATRIX  3 

DIMENSION  4(3) 

COMMON/ VAR/LlMl ,  L IM? ,  ZE,  IP,  H ,  <f  ,Dt  L 

C0MM0N/MA7/FK513)  ,F2<313)  ,F3<513),A  (3  >  *  AA  (  3  > ,  3  (s  13,  3) ,  0  (  3 ) 

A  i'  { P » Q  i  3T*3-  R  '  ' 

A2(P,Q,R> =>*R*P»P-2.*0 

A  3  TP ,  Q",  R )  si*  3 + P  *  ? '*  R  *  R  -  2 .  *  P  ♦  a  *  Q 

A*»(Rf Q*R>  s3*1 

PI=2.*3.1«159 

MH=FRlQ=0, 

oc  io  i=i  Tn  . . 

TN1={A1(H (t) »W(3) ,H(3))‘*2) 

SN=NW*WW*  HW*WW*A2(W(1)  «WC2)  «MC3)  )  +  A3(W(l)»H(2).W(3))) 

SN= (WW**c  )  +SN*  ( A*<  («J1),«(2),H(3M**2) 

TN=  (At-  ( W  (  II  *VI  (2)  *  H  (  3 ) )  *♦  2 )  *  ( WW*WW+TN1 ) 

F3( I) =TN* ?2  tl ) / (TN1#SN) 

8  (If  1)  =  AH*MM*(AA*WN*W(1)  -»W(1)  *€/»(?)♦  W<3»-2.*W (2)  I  ) 
3<Itl>sTNlM3(I»JH-A*<WCl),*M2)  »K*(3)  )*M<2)  >  *41  (W  <  1 )  ♦  W<  ? )  ♦  W  (  3  J  >*SN 
SNlsAl CMC  1) ,M(2) , W(3l)*At H(l)  .WC2) , VI 3) ) +  W 1 2) * (MW*WW+T»1 ) 
SNla(SNt*TU*A«»(M(ll«H(2)  »W(3))*SN-TN*3(I»1))+2«,F2(I> 

6 (I »1 > =SNl/<TNl*TNi *SN* jN) 

BCI.2)  =  (WH*V»M*(-HW* NW*WC2)-N(1)*W(1)  >  +  A«f<W<l>,W<2)1W(3>)*H(l))*TN 
3d,2)=(4C(N(i  )7w(2)  ,V(3)  )  *W  (1 J  ♦  <  WW*WW+Thl  J  *SM-B  < 1 ,2  ) )  *2.  *F2  ( I) 
3<I,2)=9(I,2J/(TN1*SN*SN) 

8 (I *3 )  =A1  ( <(1 )  t M J  2)  *M<3) )  ( W(3 )*WH*WM+’W(1  )*H(1I*H(3)  >-SN 

8  C I » 3 )  *-  ( A  *  (W  ( 1 1  »H  ( 2_)  iHI3l  l  **2) *TN1*SN-TN*B (l » 3  ) 

6(i»3)=2.  ♦- 2  ( I  >*  B  (I  ,3>/<  CA1(«(1»  ,W(2>,W(3)  )*»3)  *SN*SN> 

FRcQ=FR£0'H 
WW=PI»rREC  ’ 

10  CONTINUE 
>NO 


SiZ*  ?AGI  15  BSST  RUALXTT  »!£»«*** 

7950*1  jUSi/  W  AC 
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COMMENTS:  EXECUTION  OF  PROGRAM  SDCON 

1.  COMPILE  PROGRAM 

2.  INPUT:  ATTACH,  TAPE1  which  has  TRAJECTORY  INFORMATION 

•  •  ** 

®x»  ®T*  ^T’  ^T*  ^T 

3.  ATTACH  INTERNATIONAL  MATHEMATICAL  &  STATISTICAL  LIBRARY  (IMSL) 

SUBROUTINE  USED  -  VCONVO 

4 .  RUN 

5 .  OUTPUT :  TAPE2 

PARTIAL  DERIVATIVES  OF  P(2,2)  WITH  RESPECT  TO  (o^,  «2,  a3) 


PROGRAM  SCION (INPUT  »GUT3UT,TAD'il*TAPE2) 

DIMENSION  4(**09:>)  .4 1  (1020  ,3  ( 2GO ) , IWK  ( 12 >  ,C  (9, 2)  ,F2 (1024 )  ,0(1024) 
DAT  A'  ♦♦iff  ♦•1*377  ♦•1#88* 

1  4»  6  3  ♦.fci>  ,-l<*33,,«,i6fi»,l-»*  li  / 

N=102**301  L=*  04 

C  TAPE1  IS  I  l°'.‘T  TAP"  WITH  TRAJECTORY  INFO  IN  DiGPt.ES 
C "  '  T  AP£2  IS  OUTPUT  T  AP 12  WITH  PARTKL  DERIVATIVES  OF  P<2,2) 

C  IC=1  FOR  A7IMUTH  AMD  IC=2  FGR  ELEVATION 

READ*  tlC 

C  A  (I)  =0(2,  3)  Al(I)=0<2,O 

C  V=P1+P2* ( C7  DOT) * (0T  DOT ) »° 3 (OT  GCT)*(OT  DDT) 

C  F2(I)=0(I ♦I)*Q(I»I) 

DO  10  1=1 *M 
T=(I-l)*OE- 

"T(i)=-EX®  (I{3,IC>*TI*(C03(3U,IS)*T>*C(3,JC)*3IN(C(4,IC)*T>  > 
AJI)*C(1,  ID*  (EXP  (0(2,  IS)  *T)-MD) 

A1(I)=«i.XP(C(3»IC)*T)*(COS(3(t»IC)*t)*C(7,IC)*SIH(C(a,IC)*T)) 
F2JI)=A<I)*4(I)+Al(I)*Al(I)+("XP(n(2,IC)*T) )**2 

A'(I)sAl'(I)  =  F.2?,*A(I)*A(I)-C(a,I3)*A(I)*Al(I)+C(9,IC)*Al(I)*Al(I) 

REAO(l,5) CT,3(t) ,D(I),E,  :n,£OD 

IF(IC.EQ. l)GO  TO  6 

8  (I )  =EQ 

o(T$=roo 

6^ _ BJI)=F(I)  *3(1) 

D(I)=OC)  *0(1) 

10  CONTINUE 

'5  'FORMAT  (66  l’.*) 

C  PARTIAL  OF  °{?,2)  WR7  P2 

CALL  VC0NV)(A,3,N,N,IWK) 

DO  100  1=1, N 
8 (I) =0(1) 

0  (I )  =A  ( I ) 
iofl  a  (I)=ai c  > 

C  PARTIAL  0-  P (2,2)  WRT  03 

CALL  VC0NVD(A,3,N,M.:w<> 

DO  20  1=1, N 
AA=A1 (I ) 

Al( I) =A ( I ) 

"A(I)  =  AA 
20  8(I)=1. 

C  PARTIAL  OF  !>(?,?)  M’T  PI 

CALL  VCONVD(A,3,N,N.IWK> 

DO  30  1*1  ,T  ’ 

30  WRITE (2,4 1- 2(1), A  (D  ♦  'M I)  .A1  (I) 

4  '  F  CRM4T  (*»G  12  •  5 ) 

END 
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COMMENTS:  EXECUTION  OF  PROGRAM  OBS 

1.  COMPILE  ROUTINE 

2.  INPUT:  ATTACH,  TAPE2 

where  TAPE2  HAS  THE  FOLLOWING  CHANNELS: 

0T*  0T’  0T’  ^T’  ^T*  ^T  FOR  A  GIVEN  TRAJECTORY 

3.  ATTACH  D.L.  KLEINMAN  [13 ]  LIBRARY 

SUBROUTINE  USED  -  DSCRT 

4.  RUN 

5.  OUTPUT:  TIME,  MEAN  TRACKING  ERROR,  STANDARD  DEVIATION 


PROGRAM  0=3(IN®UT , OUTPUT, TA^r 2) 

DIMENSION  AU,i,),?(4|,)((,„),Rl(u),Hlli,4l,U2U,i.), 

i  P(: 5;  2)  . .  '* .  . . 

COM  MON/MA  IN  1  /  N  »  N  2 

TIME  5  0.  -  a5,  TIM-  S?c°  =  .0'* 

P  =  lAMOAI,  LAM  DA  2 ,  Ki ,  <2,_Or  DDT 
DATA  DEL  ,  N,  TEND/,  0~u  ,  E,v5  .  / 

DATA  P/-1..-1,  7.',  25.  ,2.6, 0.  ,-l.  ,- 1.3  77, 25.  ,3.76,0,./ 

N 1  =  N*  N  ~ . . .  . 

N2=N+ 1 
L 1=  1 
L2  =  2 

PPl NT ♦ , 5  0  t-  TY^f  1  R  OR  AZ'”'  2  FOR  EL  '  3  FOR  80TH 

READ*  ,IFG  _  _  _  _ 

i  f  ( i  f'g  q  .  r>  u  2=1"”  .  ' '  ‘  .  ~~  ‘  ' 

IF( IPGs  EQ.2) Li  =  2 

IC=T  FOR  A  7 1  MU  T  -T  A  N  0  ~  1 C  =T  c  0  R  ELEVATION 
00  30  0  I C  =L  1 »  L  2 
REWIND  2 

PI , P2 , °3  ARE  GOEFICIENTS  ASSOCICATED  WITH  NOISE  COVARIAN 
PRINT*,  20  H  TPf'iN  "Pi  PfP'3 
READ*, PI, P2,D3 

INITIALIZE  VECTORS  AND"  MATRICES 
00  10  1  =  1, Nl 

All)  =  0,  "  .  . 

X(I)  =  0,  _ 

m  (2)  =  l  +  P('l7TG)  . . . 

A  (5  )  =  F  (  2  ,  IS ) 

A (10) =-J <  1,10 
A(l-,)  =-0  (  2,  IS) 

A(13)=-p<  3, IS) 

A  (12)  =1. 

A  ( 1 0 )  =  - °  ( •< •  T 0 >  '  '  . . 

CALL  DSC5  T (N , A , OE  L , W1 , W 2 , 10) 

W 1= TRANS I TION  HAT=IX  W2=INTZGR£L  OF  TRANSITION  MATRIX 
JJ=1 

DO  63  J=  3  <3 
Jl=J+£ 

R 1  ( J  J  )  =  -  W  2  ( J )  -  W2  (  J 1 ) 

UJ=JJ+1 
GO  60  :=1 , J 
Z (I ) =  G*  ‘ 

DO  oO  J=1  ,  I  . 

A  (I ,  J )  =R1  ( I )  * R 1  (  J_) 

00  20  1  =  1  ",  J  ’  “  . ' 

1 1= 1+  2*N 

R1  ( I)  =  W2  (  I)  +  W2  ( I  I  > *  <  1-P  ( 5”,  IC >  ) 

T=0« 

READ  IN  TRAJECTORY  INpO  FOR  AT  AND  EL  IN  DEGREES 
REA 0(2, 3)  (P (1,1)  ,  1=1,3) , (P(I, 2) ,1  =  1,3) 

FORMAT  (6G  12,  -*) 

COMPUTE  ESrIMATEr  OT  POT  AND  OT  C DT 


3313  PASS  IS  SS-r- 


I  §UALm 


DOO 
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3  0 


c 

c 


15 

25 


35 

C 


40 


C 

C 

C 


75 

500 


5 

10 


20 


30 


P4  =Z(3)-Z  <1> 

Pf=(P—PP<,)/3EL 
1F(T.z3.0.)P5s3(?,:C) 
v=noise  ccvariavo;: 
Va(3i*oU4p+*32*P5*P5+P3l/0iu 
PPv*P«, 

COMPUTE  MEAN  STATE  EQUATION 

00  25  1=1  ,M 

W2(I>=0. 

11*1 

DO  Is  J=I  » 'll » N 
W2< II =W2 { II+H1IJI 

ii*r:n 

CONTINJI 
DO  35  1=1 , I 

Z  (I ) =W2 ( I  )»Rl (I) ♦ P( 3 , ICI 

CONTINUE  ‘  - -  - - 

COMPUTE  COVARIA  1C  *  MATRIX  JF  STA'u  EQUATION 
CALL  MULT  Hi, X.‘I,N1,H?)  ' 

oc  4j  i*i, n 

X  (I)=A(I)  *v+w?('i)  *  ' 

SG=SQRT {X <’,?)) 

T*T  +OtL  . 

MEAN  TRACKING  ERROR  IS  7(21 
VARIANCE  C-  TRACKING  fP.RO0  IS  X  <2,21 
PRINT  OUTFJT  AT  EVERT  SECOND 
LK=(T-r.0Gl>/0rL 

IF(M3C(L<  ,25) .  £0.  0)  PRINCE  ,T, Z(2) ,S0 

FORMAT (3X,TG12;,)  . . 

iF(T.GE.T-IO)GO  TO  500 

GO  TO  1 

CONTINUE 

END 

SUBROUTINE  MJLTti ,r ,L,L1,H) 

DIMENSION  EtLll.'p  <Lf)VG(16f;H(U) 

DO  10  1  =  1, L 
11=1 

00  13  <= 1 »— 

TEMP=3. 

DO  5  J= jl  ,  LI ,  L 

TSM?=TEMP+E<J>‘*tm  '  "  '  ~ 

1X»XX*1 
KK=«-1)*L*I 
G(<<)  =T“MF 

DO  23  2  =  1  ,.  '  "  -- 

00  20  S*I,fc 

TLMP=0,'  '  * . 

1I=K 

DO  15  J=I  ,ul ,u  5 

TEMP* TIM® ♦G (J)*E(II) 

11= II+L 
KK=«-1)*L*I 

H IKK) =T  £M  P  . 

L2=L-1 

UC  33  1=1, l2 
L3=I*1 

00  30  J=L  3, L 
K1=<I-1)*H-J 
K2=(J-U*UI  '  ' 

H(K1)=MI<2) 

cNO  -  -  - 


yj? 

A  & 

& 


i 
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